Tutorial 9: Advanced Modelling and Visualisation

Programming for Engineers (ENGG1001)

Tutorial Tasks

• Task 1: Creating a Meshgrid for 3D Surface Plots
• Task 2: Creating a Planar Surface
• Task 3: Modelling waves with a Radial Basis Function
• Task 4: Visualising 3D Surfaces
• Task 5: Identifying local maxima and minima
• Task 6: Animating a 3D Wave Surface

Task 1: Creating a Meshgrid for 3D Surface Plots

So far, you have had practice with creating 1D numpy arrays and plotting them. For 3D surface plots, we need to create a grid of points in the x-y plane. Write a function create_meshgrid(length: float, num_points: int) -> tuple[np.ndarray, np.ndarray] that creates a meshgrid of points from -length to length with num_points for both x and y axes. For a further explanation of what a meshgrid is, see the callout below.

NoteDefinition

What exactly is a meshgrid? A meshgrid is a 2D grid of points that can be used for plotting surfaces or evaluating functions over a 2D domain. It consists of two 2D arrays. Consider the example below where we create a 10 x 10 meshgrid:

x_mesh:
[[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]]
y_mesh:
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
 [3. 3. 3. 3. 3. 3. 3. 3. 3. 3.]
 [4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]
 [5. 5. 5. 5. 5. 5. 5. 5. 5. 5.]
 [6. 6. 6. 6. 6. 6. 6. 6. 6. 6.]
 [7. 7. 7. 7. 7. 7. 7. 7. 7. 7.]
 [8. 8. 8. 8. 8. 8. 8. 8. 8. 8.]
 [9. 9. 9. 9. 9. 9. 9. 9. 9. 9.]]

Essentailly, a meshgrid is a way of presenting a set of cartesian coordinates in a structured format. So if we feed an index (which essentially acts as a coordinate) into the meshgrid, we can get the corresponding x and y values. For example, if we want to get the x and y values at the index (2, 3), we can do:

x value at x_mesh[2, 3]: 3.0
y value at y_mesh[2, 3]: 2.0

Here is a sample usage of the function:

>>>

x_mesh, y_mesh = create_meshgrid(10, 100)

>>>

print(x_mesh.shape)

 

(100, 100)

>>>

print(x_mesh[0, 0])

 

-10.0

>>>

print(y_mesh[50, 0])

 

3.303030303030301

Tip

What function could you use to create the 1D arrays for x and y? Once you have created the 1D arrays, then you can use numpy.meshgrid to create the 2D grid.

Callout-answerAnswer

meshgrid.py

import numpy as np

def create_meshgrid(length: float, num_points: int) -> np.ndarray:
    """Creates a meshgrid of points from -length to length with num_points for both x and y axes.

    Parameters
    ----------
    length : float
        The half-length of the grid in both x and y directions.
    num_points : int
        The number of points to generate along each axis.
    """
    x = np.linspace(-length, length, num_points)
    y = np.linspace(-length, length, num_points)
    x_mesh, y_mesh = np.meshgrid(x, y)
    return x_mesh, y_mesh

x_mesh, y_mesh = create_meshgrid(10, 100)
print(x_mesh.shape)
print(x_mesh[0, 0])
print(y_mesh[50, 0])

Output

(100, 100)
-10.0
0.10101010101010033

Task 2: Creating a Plane Surface

Now that we can create a meshgrid of points, we can use it to create a 3D surface. Write a function plane_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, a: float, b: float, c: float) -> np.ndarray that computes the z values for a plane surface defined by the equation: \[z = ax + by + c\]

Here is a sample usage of the function:

>>>

x_mesh, y_mesh = create_meshgrid(10, 100)

>>>

plane = plane_surface(x_mesh, y_mesh, 1.0, 2.0, 3.0)

>>>

print(plane.shape)

 

(100, 100)

>>>

print(plane[0, 0])

 

-27.0

>>>

print(plane[50, 50])

 

3.303030303030301

Tip

To compute the wave surface, you will need to calculate the radial distance R for each point in the meshgrid and then apply the RBF formula. Remember to use vectorisation!

Callout-answerAnswer

plane_surface.py

def plane_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, a: float, b: float, c: float) -> np.ndarray:
    """Computes the z values for a plane surface defined by the equation z = ax + by + c.

    Parameters
    ----------
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    a : float
        The coefficient for x.
    b : float
        The coefficient for y.
    c : float
        The constant term.

    Returns
    -------
    np.ndarray
        The computed z values for the plane surface.
    """
    z = a * x_mesh + b * y_mesh + c
    return z
x_mesh, y_mesh = create_meshgrid(10, 100)
plane = plane_surface(x_mesh, y_mesh, 1.0, 2.0, 3.0)
print(plane.shape)
print(plane[0, 0])
print(plane[50, 50])

Output

(100, 100)
-27.0
3.303030303030301

Task 3: Modelling waves with a Radial Basis Function

Now that we have a meshgrid of points, we can model a wave surface using a Radial Basis Function (RBF). The RBF is defined as: \[f(x, y) = A \cdot \cos(kr - \omega t) \cdot e^{-\alpha r}\] Where:

• \(A\) is the amplitude of the wave
• \(k\) is the wave number
• \(r\) is the radial distance from the origin: \(r = \sqrt{x^2 + y^2}\)
• \(\omega\) is the angular frequency
• \(\alpha\) is the decay constant of the wave with respect to distance
• \(t\) is time in seconds

Write a function radial_wave(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float, k: float, omega: float, alpha: float, t: float) -> np.ndarray that computes the wave surface values for the given parameters. Here is a sample usage of the function:

>>>

x_mesh, y_mesh = create_meshgrid(10, 100)

>>>

wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)

>>>

print(wave_surface.shape)

 

(100, 100)

>>>

print(wave_surface[50, 50])

 

0.9458561801297212

>>>

print(wave_surface[0, 0])

 

-0.24310473049518086

Callout-answerAnswer

radial_wave.py

def radial_wave(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float, k: float, omega: float, alpha: float, t: float) -> np.ndarray:
    """Computes the wave surface values for the given parameters using a Radial Basis Function.

    Parameters
    ----------
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    A : float
        The amplitude of the wave.
    k : float
        The wave number.
    omega : float
        The angular frequency.
    alpha : float
        The decay constant of the wave with respect to distance.
    t : float
        Time in seconds.

    Returns
    -------
    np.ndarray
        The computed wave surface values.
    """
    r = np.sqrt(x_mesh**2 + y_mesh**2)
    wave_surface = A * np.cos(k * r - omega * t) * np.exp(-alpha * r)
    return wave_surface
x_mesh, y_mesh = create_meshgrid(10, 100)
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
print(wave_surface.shape)
print(wave_surface[50, 50])
print(wave_surface[0, 0])

Output

(100, 100)
0.9458561801297212
-0.24310473049518086

Task 4: Visualising 3D Surfaces

Now that we have a function to compute the wave surface, we can visualise it using Matplotlib’s 3D plotting capabilities. Write a function plot_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, surface: np.ndarray) that creates a 3D surface plot of the surface.

Here is a sample usage of the function potting the plane surface:

>>>

x_mesh, y_mesh = create_meshgrid(10, 100)

>>>

plane = plane_surface(x_mesh, y_mesh, 1.0, 2.0, 3.0)

>>>

plot_surface(x_mesh, y_mesh, plane)

Which produces the following 3D surface plot:

3D Surface Plot of the Plane Surface

Also try plotting the wave surface using the same function:

>>>

x_mesh, y_mesh = create_meshgrid(10, 100)

>>>

wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)

>>>

plot_surface(x_mesh, y_mesh, wave_surface)

Which produces the following 3D surface plot:

3D Surface Plot of the Radial Wave

Tip

To be able to setup an axis for a 3D plot, you need to specify the projection argument when creating the subplot. The plot_surface method can be used to create the surface plot. You can use the code below to get started:

fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
ax.plot_surface(X, Y, Z)

If you want to make your plot look nicer, you can also specify a colormap using the cmap argument in plot_surface. For example, you can use cmap='viridis' to use the Viridis colormap, like the example above.

Callout-answerAnswer

plot_surface.py

import matplotlib.pyplot as plt

def plot_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, surface: np.ndarray):
    """Creates a 3D surface plot of the surface.

    Parameters
    ----------
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    surface : np.ndarray
        The computed surface values.
    """
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
    ax.plot_surface(x_mesh, y_mesh, surface, cmap='viridis')
    ax.set_xlabel('X-axis')
    ax.set_ylabel('Y-axis')
    ax.set_zlabel('Surface')
    ax.set_title('3D Surface Plot')
    plt.show()

x_mesh, y_mesh = create_meshgrid(10, 100)
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
plot_surface(x_mesh, y_mesh, wave_surface)

Output

Note

Now that you have created an easy way to visualise the wave surface, determine what effects the following parameters have on the wave surface by changing their values and observing the resulting plot:

1. Amplitude A
2. Wave number k
3. Angular frequency omega
4. Decay constant alpha
5. Num points in the meshgrid num_points

Task 5: identifying local maxima and minima

Now that we have a visualisation of the wave surface, we can identify the local maxima and minima of the surface. Write a function find_local_extrema(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray] that finds and returns the radius values (how far from the origin) of the local maxima and minima of the surface. Here is a sample usage of the function:

>>>

x_mesh, y_mesh = create_meshgrid(10, 100)

>>>

wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)

>>>

maxima, minima = find_local_extrema(wave_surface)

>>>

print(maxima)

 

[9.39448244 6.16244406 3.1329419 0.14284985 3.1329419 6.16244406

 

9.39448244]

>>>

print(minima)

 

[7.77843366 4.7485492 1.51851479 1.51851479 4.7485492 7.77843366]

Tip

To find the local maxima and minima of radial function, we can use the fact that the surface is “radially symmetric”, meaning that the surface values only depend on the radius from the origin. As such, we can compute the radius values for a central “slice” of the surface (e.g. a given row or column of the surface), and then find the local maxima and minima of that slice.

To illustrate this, consider the following central slice of the curve above:

Central Slice of the Radial Wave Surface

We can see that the local maxima and minima of the surface correspond to the local maxima and minima of this slice. Therefore, we can find the local maxima and minima of this slice to determine the radius values of the local maxima and minima of the surface.

To find the local maxima and minima of a 1D array, you can start by using the np.diff function to compute the differences between consecutive elements. Then you can find where the sign of these differences changes (i.e. from positive to negative for maxima, and from negative to positive for minima). Then you can store the corresponding radius values for those indices and return them as the output of the function.

Callout-answerAnswer

Two solutions are presented below, one using a more manual approach, and the other using vectorisation:

find_local_extrema.py

def find_local_extrema_simple(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Finds and returns the radius values of the local maxima and minima of the surface.

    Parameters
    ----------
    surface : np.ndarray
        The computed surface values.

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        A tuple containing two arrays: the first array contains the radius values of the local maxima, and the second array contains the radius values of the local minima.
    """
    # Compute the radius values for a slice of the 
    # surface (e.g. the middle row)
    radius = np.sqrt(x_mesh[x_mesh.shape[0] // 2]**2 
                    + y_mesh[y_mesh.shape[0] // 2]**2)
    
    diffs = np.diff(surface[surface.shape[0] // 2])  # taking the middle row as the slice

    maxima = []
    minima = []
    for i in range(1, len(diffs) - 1):
        if diffs[i-1] > 0 and diffs[i] <= 0:  # local maximum
            maxima.append(radius[i])
        elif diffs[i-1] < 0 and diffs[i] >= 0:  # local minimum
            minima.append(radius[i])
    
    return np.array(maxima), np.array(minima)

def find_local_extrema(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Finds and returns the radius values of the local maxima and minima of the surface.

    Parameters
    ----------
    surface : np.ndarray
        The computed surface values.

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        A tuple containing two arrays: the first array contains the radius values of the local maxima, and the second array contains the radius values of the local minima.
    """
    # Compute the radius values for a slice of the 
    # surface (e.g. the middle row)
    radius = np.sqrt(x_mesh[x_mesh.shape[0] // 2]**2 
                    + y_mesh[y_mesh.shape[0] // 2]**2)
    
    # Compute the differences between consecutive elements in the surface slice
    diff = np.diff(surface[surface.shape[0] // 2])  # taking the middle row as the slice
    
    # Find where the sign of the differences changes to identify local maxima and minima
    maxima_indices = np.where((diff[:-1] > 0) & (diff[1:] <= 0))[0] + 1
    minima_indices = np.where((diff[:-1] < 0) & (diff[1:] >= 0))[0] + 1
    
    # Get the corresponding radius values for those indices
    maxima = radius[maxima_indices]
    minima = radius[minima_indices]
    
    return maxima, minima
x_mesh, y_mesh = create_meshgrid(10, 100)
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
maxima, minima = find_local_extrema(wave_surface)
print(maxima)
print(minima)

Output

[9.39448244 6.16244406 3.1329419  0.14284985 3.1329419  6.16244406
 9.39448244]
[7.77843366 4.7485492  1.51851479 1.51851479 4.7485492  7.77843366]

Task 6: Animating a 3D Wave Surface

For a challenge task, try to create an animation of the wave surface over time. However there is a slight twist: the equation for the wave surface now accounts for time decay as well, so the wave surface will not only change with time, but it will also decay over time. The new equation is: \[f(x, y, t) = A \cdot \cos(kr - \omega t) \cdot e^{-\alpha r} \cdot e^{-\beta t}\] Where \(\beta\) is the decay constant of the wave with respect to time. You can use Matplotlib’s animation capabilities to create the animation. Here is a sample usage of the function:

>>>

x_mesh, y_mesh = create_meshgrid(10, 100)

>>>

animate_wave_surface(x_mesh, y_mesh, A=1.0, k=2.0,

...

omega=1.0, alpha=0.1, beta=0.05,

...

duration=10, fps=30)

Figure 1: Animated 3D Wave Surface Plot

Tip

Use the following incomplete code as a starting point for your animation function. You will need to complete the update_image function and the animate_wave_surface function to create the animation.

from matplotlib import animation
def update_image(frame, x_mesh, y_mesh, A, k, omega, alpha, beta, ax):
    """Updates the wave surface for each frame of the animation.

    Parameters
    ----------
    frame : int
        The current frame number.
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    A : float
        The amplitude of the wave.
    k : float
        The wave number.
    omega : float
        The angular frequency.
    alpha : float
        The decay constant of the wave with respect to distance.
    beta : float
        The decay constant of the wave with respect to time.
    ax : matplotlib.axes._subplots.Axes3DSubplot
        The 3D axis object to update.
    """
    wave_surface = # write your code here!

    ax.clear()
    ax.plot_surface(x_mesh, y_mesh, wave_surface)
    ax.set_xlabel('X-axis')
    ax.set_ylabel('Y-axis')
    ax.set_zlabel('Wave Surface')
    ax.set_title('Animated 3D Wave Surface Plot')

    return [ax]



def animate_wave_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float,
                         k: float, omega: float, alpha: float, beta: float,
                         delta_t: float, num_steps: int):
    """Creates an animation of the wave surface over time.

    Parameters
    ----------
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    A : float
        The amplitude of the wave.
    k : float
        The wave number.
    omega : float
        The angular frequency.
    alpha : float
        The decay constant of the wave with respect to distance.
    beta : float
        The decay constant of the wave with respect to time.
    delta_t : float
        The time step for each frame of the animation.
    num_steps : int
        The number of steps in the animation.
    """
    
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

    #| decorate the plot with labels and title here!

    anim = animation.FuncAnimation(fig, update_image, frames=np.arange(0, num_steps, delta_t),
                                   fargs=(x_mesh, y_mesh, A, k, omega, alpha, beta, ax), interval=50)
    plt.show()

Callout-answerAnswer

animate_wave_surface.py

from matplotlib import animation

def update_image(frame, x_mesh, y_mesh, A, k, omega, alpha, beta, ax):
    """Updates the wave surface for each frame of the animation.

    Parameters
    ----------
    frame : int
        The current frame number.
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    A : float
        The amplitude of the wave.
    k : float
        The wave number.
    omega : float
        The angular frequency.
    alpha : float
        The decay constant of the wave with respect to distance.
    beta : float
        The decay constant of the wave with respect to time.
    ax : matplotlib.axes._subplots.Axes3DSubplot
        The 3D axis object to update.
    """
    R = np.sqrt(x_mesh**2 + y_mesh**2)
    wave_surface = A * np.cos(k * R - omega * frame) * np.exp(-alpha * R) * np.exp(-beta * frame)

    ax.clear()
    ax.plot_surface(x_mesh, y_mesh, wave_surface)
    ax.set_xlabel('X-axis')
    ax.set_ylabel('Y-axis')
    ax.set_zlabel('Wave Surface')
    ax.set_title('Animated 3D Wave Surface Plot')

    return [ax]

def animate_wave_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float,
                         k: float, omega: float, alpha: float, beta: float,
                         delta_t: float, num_steps: int):
    """Creates an animation of the wave surface over time.

    Parameters
    ----------
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    A : float
        The amplitude of the wave.
    k : float
        The wave number.
    omega : float
        The angular frequency.
    alpha : float
        The decay constant of the wave with respect to distance.
    beta : float
        The decay constant of the wave with respect to time.
    delta_t : float
        The time step for each frame of the animation.
    num_steps : int
        The number of steps in the animation.
    """
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
    ax.set_xlabel('X-axis')
    ax.set_ylabel('Y-axis')
    ax.set_zlabel('Wave Surface')
    ax.set_title('Animated 3D Wave Surface Plot')

    anim = animation.FuncAnimation(fig, update_image, frames=np.arange(0, num_steps, delta_t),
                                   fargs=(x_mesh, y_mesh, A, k, omega, alpha, beta, ax), interval=50)
    plt.show()